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We report a characterization of the relative stability and structural behavior of various micellar 
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crystals of an athermal model of AB-diblock copolymers in solution. We adopt a previously devel- 
oped coarse-graining representation of the chains which maps each copolymer on a soft dumbbell. 
Thanks to this strong reduction of degrees of freedom, we are able to investigate large aggregated 
systems, and for a specific length ratio of the blocks / = Ma/ (Ma + Mg) = 0.6, to locate the 

O 

order-disorder transition of the system of micelles. Above the transition, mechanical and thermal 
properties are found to depend on the number of particles per lattice site in the simulation box, 
and the application of a recent methodology for multiple occupancy crystals (B.M. Mladek et al., 
Q_! Phys. Rev. Lett. 99, 235702 (2007)) is necessary to correctly define the equilibrium state. Within 

this scheme we have performed free energy calculations at two reduced density p/p* = 4, 5 and 
for several cubic structures as FCC,BCC,A15. At both densities, the BCC symmetry is found to 

O 

correspond to the minimum of the unconstrained free energy, that is to the stable symmetry among 
the few considered, while the A15 structure is almost degenerate, indicating that the present sys- 
tern prefers to crystallize in less packed structures. At p/p* = 4 close to melting, the Lindemann 
ratio is fairly high (~ 0.29) and the concentration of vacancies is roughly 6%. At p/p* = 5 the 
mechanical stability of the stable BCC structure increases and the concentration of vacancies ac- 
cordingly decreases. The ratio of the corona layer thickness to the core radius is found to be in 
good agreement with experimental data for poly(styrene-b-isoprene) (22-12) in isoprene selective 
solvent which is also reported to crystallize in the BCC structure. 
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I. INTRODUCTION 



A recurrent phenomenology in soft matter is the self-assembly of its basic constituents into 
supramolecular aggregates of various shapes and morphologies. This property, ubiquitous 
in materials of biological and technological interest, greatly enriches the phase behavior of 
these systems with the possibility of order-disorder transition of aggregates on mesoscopic 
time and space scales. A remarkable example in this sense is represented by diblock copoly- 
mers both in melts [H [2] and solution [3j Ej. The occurrence of micro-phase separations into 
a large number of spatially patterned arrangements (lamellar, gyroid) or clusters of different 
shape and size (cylinders, micelles) in diblock copolymer solutions has been experimentally 
characterized [3HH]- For clustered conditions, it has also been observed the occurrence of a 
transition between a disordered phase of dispersed clusters and a phase with long range spa- 
tial order in hexagonal (for cylindrical aggregates) or cubic (for spherical clusters) structures 
(ODT). Moreover thermotropic and lyotropic transitions between cubic crystals of different 
structure have been reported [7HTU] and theoretically analyzed by mean field theoretical 
arguments [TTHT3] . 

Thanks to the number of available control parameters, such as temperature, concentration 
of the copolymer, solvent selectivity with respect to the two blocks and the relative block 
lengths, these systems are in principle very versatile and could be designed to obtain a specific 
phase behavior. In practice however this program requires a microscopic understanding of 
the entropic-enthalpic competition behind the self-assembling and the ODT which is still 
missing. This complex picture requires a considerable effort both in theoretical and molecular 
simulation approaches, and while in the melt regime the low osmotic compressibility allows 
to successfully apply mean field theory [TJ [2] , in solutions the possibility to make ab initio 
predictions is more difficult pj)] . Even for implicit solvent models the application of the 
common molecular simulation strategies is limited by the huge number of degrees of freedom 
in the system (many long chains) and the wide range of time and length scales involved [T5l 
[16]. In order to describe the phase behavior of large aggregated systems it is tempting to 
use a coarse graining scheme which allows to reduce the number of degrees of freedom by 
mapping a number of monomers of each chain onto soft blobs. This strategy goes back to 
a seminal work of Louis et al. [17] which mapped a solution of homopolymers onto a fluid 
of soft colloids interacting through, density dependent, soft core, pair potentials obtained 
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from the observed structure through liquid state theory methods. Working with density 
dependent potentials is however cumbersome 1 1 8] . in particular for inhomogeneous systems, 
and can be avoided by switching to the so called "multi-blob" representation, in which the 
degree of coarse graining (the number of monomer mapped onto a single blob) is tuned on 
the local density. This scheme has been applied to a semi-dilute solution of homopolymers 
in good solvent [19] and very recently an extension to diblock copolymer solutions has been 
proposed[20j. 

Within this strategy, the coarsest grained representation of diblock copolymers maps 
each block of a chain onto its center of mass reducing a long polymer to a soft dumbbell. 
This is the minimal model of diblock copolymers and it is accurate at high dilution only. 
Its extension to finite concentration (i.e. determining density dependent potentials as in 
the case of homopolymers) is not trivial since the theoretical connection between structure 
and pair interactions at finite density is missing. Nonetheless the minimal model with zero 
density interactions, extensively investigated by Monte Carlo (MC) simulations well beyond 
the overlap concentration [2TH23] , exhibits the self-assembling phenomena and a number of 
ODT transitions in qualitative agreement with the experimental phenomenology of diblock 
copolymer solutions. This model has been studied for various solvent selectivity and over a 
wide range of densities and block lengths ratios [2~Ij 123"]. It has been found that the ISI model, 
in which each block of the copolymer is in an ideal (theta) solvent while pairs of monomers 
of the two blocks repel each other (mutually-avoiding), has a microphase separation towards 
a lamellar phase beyond a given concentration [21], whereas a ISS model, in which one 
block is in theta solvent while the other is in a good solvent and pairs of monomers of 
differents block are mutually-avoiding, naturally form spherical micelles beyond a CMC 
density [22], [23] . Moreover in the ISS model, an order-disorder transition (ODT) between a 
liquid and a crystalline phase of micelles is spontaneously observed to occur for high enough 
concentration. The same qualitative picture emerges from a "multi-blob" investigation of 
the ISS model [20] with the addition of an observed cylindrical phase. Because of the large 
aggregation number (~ 100), the number of micelles in a simulation sample is limited and 
the structure of the crystals that spontaneously form is influenced by the constraints of the 
simulation and by the slow relaxation of the system. In order to obtain detailed predictions 
about the relative stability of disordered and ordered phases of different symmetries, free 
energy calculations are necessary. 
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The computational approach to the free energy of these complex systems is significantly 
more difficult than the usual methods employed in simple non-assembling systems. On 
the one hand, the soft interactions between blobs makes the use of Widom insertion method 
highly efficient even at high density and even in the crystal phase (the phase in which micelles 
are ordered according to a lattice). On the other hand, the softness of the potential produces 
floppy and polydisperse aggregates which pose a fundamental problem to the calculation of 
the free energy of crystalline phases. The problem is related to how correctly to define 
the equilibrium state of the system with broken symmetry, in particular to determine the 
equilibrium number of particles per lattice site for a given density (and temperature in the 
general case). This problem does not arise in crystals with single-particle occupancy (one 
particle per lattice site) but has been shown to be relevant for multiple-particle occupancy 
crystals [SJ [25]. For the minimal dumbbell model of dibock copolymers beyond the ODT, 
the size distribution of the micelles is found to depend significantly on the number of particles 
in the sample and a given crystal structure with fixed number of lattice sites can comfortably 
accommodate a wide range of particle numbers, by adjusting the size of the aggregates, 
before loosing its structure. In a recent paper on multiple-occupancy crystals, using ideas 
from the thermodynamic theory of crystals with vacancies [2B], Mladek et al. [2~4"l |2"5] devised 
a method to define the equilibrium number of particles per lattice site and a strategy to 
compute the free energy of the equilibrium system. In the present paper we apply the same 
formalism and method to compute the free energy of cubic crystals of spherical micelles of 
diblock copolymers in solution. Since the method is quite demanding in terms of computer 
time, we adopt the dumbbells minimal model of diblock copolymers despite the fact that 
its use at finite density is not fully justified. Future work will consider the more accurate 
multiblob representation. 

Finally, we should mention a previous attempt to address the relative stability of crystals 
of different symmetries of diblock copolymers micelles [27]. In this approach micelles are 
mapped onto point particles interacting through density dependent effective pair potentials 
obtained by inverting the micelle-micelle radial distribution function. This work suggests 
that in the range 0.4 < / < 0.6 the OD transition density decreases upon increasing /, and 
that the A15 structure [28j, already proposed as the stable structure for spherical micelles 
in linear block copolymer melts [311 E2] , in start block copolymer melts [33j H2] and exper- 
imentally observed for dendrimic clusters [35], is the most stable structure among the few 
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explored, namely BCC, FCC, diamond, SC, A15. This is somehow in contradiction with 
the expectation that the favored structure should depend on f — Ma/ (Ma + Mb), with 
micelles with smaller coronae (crew-cut) preferring close packed structures such as FCC 
or HCP, while micelles with larger coronae (hairy) prefer less packed structures such as 
BCC and A15 [5J El H2]- However this effective micelle model differs from the original sys- 
tem of dumbbell micelles in several important aspects. At the dumbbell level, the micelles 
are highly polydisperse presenting a gaussian-like aggregation number distribution peaked 
around 100 dumbbells with a width of roughly 40 dumbbells[22, 23J. The polydispersity is 
lost in mapping micelles onto point particles. Moreover, the number of aggregates in the 
system at the dumbbell level of description is found to fluctuate producing a crystal with 
a finite concentration of vacancies, this might be mainly due to the floppy character of the 
aggregates. Finally, even at densities beyond the observed ODT, in the system of dumb- 
bells some molecules remain dispersed in the volume between the clusters, establishing a 
"chemical" equilibrium between the aggregated and non aggregated molecules. Again this 
phenomenon is lost when mapping micelles onto point particles, although it is implicitly 
taken into account in the effective interactions. Because of such important differences the 
prediction of the "effective" micelles system can be only considered as qualitative and a 
definitive assessment should be provided by a direct calculation of crystal free energies at 
the dumbbell level of description. 

The paper is organized as follows. In section [IT] we briefly review our model which has 

we 



been already presented in previous works [22, 23J. In the following sections III and IV 



review the thermodynamic theory of cluster crystals and the method to compute the crystal 
free energy. We also describe specific aspects related to our present system. In section 
fy] we report results for the disordered phase and the approach to the spontaneous OD 
transition which were not completely discussed in previous publications. In the following 
section we discuss our results of the free energy calculations of several crystalline structures 



at two different densities. Finally in section VII we draw our conclusions and perspectives. 
Appendices I and II are devoted to the calculations of the free energy of the reference system, 
and to details about the application of the Gauss-Legendre method to perform the coupling 
constant integration, respectively. 
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II. COARSE-GRAINED MODEL 



Let us consider a system of N copolymer chains, each formed by two blocks A and B of 
Ma and Mb monomers of different chemical species, in a volume V at fixed chain density 
p = N/V. We consider an implicit solvent model in which solvent degrees of freedom are 
absent and the solvent selectivity is represented by effective monomer-monomer interactions 
specific for each different pair of particles. We consider here the ISS athermal model repre- 
senting theta solvent conditions for monomers of the A block and good solvent conditions for 
monomer of the block B. Moreover the interaction between monomers of A and B blocks is 
mutually-avoiding. The phase diagram of this model is fully determined by entropic effects 
and the only independent variables are the chain density p and the relative size of the blocks 
/ = Ma/ {Ma + M B ). As usual in polymer solutions [29], the density is better expressed in 
terms of the reduced density p/p* where p* = 3/(4TrRg) is the overlapping density and R g is 
the radius of gyration of an isolated chain. Direct simulation of this full monomer model in 
the scaling regime (i.e. for long enough chains) is challenging because of the large number of 
degrees of freedom in the system (N x M) and the wide range of relaxation times involved. 

The simplest coarse grained representation of this model is obtained by mapping each 
block on its centre of mass therefore reducing each chain to a simple dumbbell [2TH25] . The 
interaction potentials between the two blobs of the same dumbbell 7ab(^) and between blobs 
of different dumbbells, vaa{t), vab(j), vbb{t) can in principle be obtained by inverting the 
structure obtained in a full monomer simulation. This is a well known procedure for single 
and multicomponent systems of homopolymers [30] . For diblock copolymers the connectiv- 
ity constraint makes the problem well defined only at zero density, i.e. for two isolated 
dumbbells [2Tj. In this case the block center of mass radial distribution functions, obtained 
by Monte Carlo simulation of a system of two interacting chains at the full monomer level, 
can be inverted to get the pair interaction potentials between the blobs [21]. However it is 
not known how to consistently extend this procedure to finite density. Therefore to study 
the interesting thermodynamic behavior of this system at finite density we must rely on the 
approximation of using density independent potentials identical to the ones extracted at zero 
density. As already observed [22J, these potentials do not differ much from the potentials 
obtained for mixtures of untethered chains at zero density. In the absence of the tethering 
constraint between A and B blocks we have a mixture of two chemically different homopoly- 



6 



mers in the presence of solvent of different quality. In this case the multicomponent HNC 
closure could be exploited to invert the pair structure and therefore extract pair potentials 
between blobs (center of mass of the homopolymers) even at finite density. However the 
physical behavior of this mixture at finite density is radically different from the behavior 
of the diblock copolymer solutions. Indeed in the mixture a macroscopic phase separation 
occurs at large enough density while in the copolymer system only microphase separations, 
such as micellization, occurs. It is therefore clear that the analogy between the two sets of 
potentials cannot be extended to finite density. 

As in previous investigations [221 123] we adopt here the potentials for a mixture of A and 
B homopolymers at zero concentration. We assume that the A chains are ideal (theta solvent) 
and the B chains are self-avoding (good solvent) and that the interaction between monomers 
of the A and the B chains are mutually-avoiding. At zero density the potentials between 
the centres of mass of any pair of such polymers v a p(r), (a, (3 = A,B) are the potentials of 
mean force 

v a p{r) = -k B T\ogg a p(r) (1) 

where g a /3{r), the centre of mass radial distribution functions, are obtained by MC simula- 
tions of a systems of two chains on lattice at the full monomer level. Morever, in order to 
restore the connectivity constraint in the diblock copolymer system, the tethering potential 
7ab(7") between the A and B blocks of the same copolymer, is consistently obtained at zero 
density from the radial distribution function between the centre of mass of the two blocks 
sab(t) obtained by MC simulation of a single diblock chain at the full monomer level 

7ab(0 = -k B T\ogs A B{r) (2) 

The potentials are reported in Fig. [TJ The effective intermolecular potential VAA( r ) — 
while vab(t) and fBsO") are s °ft and roughly gaussian with a characteristic length compa- 
rable to the copolymers gyration radius, and an amplitude of roughly 2kbT at full overlap. 
We make the assumption that the zero density potentials can be used at finite density. This 
transferability assumption is supported by the evidence of a weak density dependence of the 
effective potentials for homopolymers [T7] . 

In previous works [2"2l |2"3] the thermodynamic behavior of this model has been explored. 
Beyond some value of the reduced density which depends on /, the dumbbells self- assemble 
into spherical micelles. At the same time thermodynamic properties such as the excess pres- 
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FIG. 1: f=0.6: effective pair potentials vaa, vab, vbb,^ab, as a function of the distance expressed 
in gyration radius units R g . 

sure and the excess chemical potential exhibit a change of behavior providing a macroscopic 
signature of the CMC and of the microphase separation. Upon further increase of the re- 
duced density a ODT is observed with a spontaneous breaking of the translation symmetry 
and an ordering of micelles according to some crystalline structure. The precise location of 
the CMC and of the ODT is obviously influenced by the zero density potential approxima- 
tion, but the qualitative picture is confirmed by the recent multi-blob exploration [20]. 

To complete the thermodynamic characterization of the dumbbell model it would be de- 
sirable to acquire knowledge about the stable crystalline structures that one should expect 
beyond the ODT. Indeed different crystalline phases under different thermodynamic con- 
ditions have been reported in experiments [3j and a number of experimental studies about 
solid-liquid and solid-solid phase transition in diblock-copolymer micelles have appeared 
[31 El El E] ■ It is tempting to explore the capability of our simple model to qualitatively re- 
produce such transitions. As mentioned above however free energy calculations for micellar 
crystals present some inherent difficulty related to the definition of the true equilibrium state 
in a simulation system with Periodic Boundary Conditions. Naively one could think that the 
free energy of a system of micelles for a given / and at given reduced density p/p*, ordered 
along the nodes of a specific crystalline lattice, can be simply obtained by computing the 
single molecule chemical potential by the Widom insertion method, which remains very effi- 
cient even at relatively high density, and the pressure by the virial expression. However our 
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floppy micelles are able to accommodate a variable number of dumbbells and in a periodic 
system at fixed density and fixed number of micelles (each one at a different lattice node) 
the number of molecules and the corresponding volume can vary in a finite range before the 
crystalline state becomes dynamically unstable. This anomaly is a clear signature of the 
presence of an additional hidden variable. As explained in the next section this variable is 
the number of lattice nodes which, in a system with Periodic Boundary Conditions (PBC), 
plays the role of a thermodynamic constraint. An analogous situation has been recently 
found in the multiple occupancy crystals, i.e. systems which in the crystalline phase are 
able to accommodate an increasing number of particles per lattice site upon increasing the 
external pressure [2U [25] . 



III. THERMODYNAMICS OF CLUSTER CRYSTALS 

Our presentation of the theory and of the strategy to compute the crystal free energies 
follows closely the original work of Swope and Anderson [2B] and a recent implementation 
by Mladek et al. j2U [25]. In the study of single occupancy crystals, it is common practice 
to associate the equilibrium state of the system with the defect-free state, except when 
one is interested in the properties of the defects. Indeed in most cases the equilibrium 
concentration of vacancies or defects is so low that this association is quantitatively correct 
[55] . For multiple occupancy or cluster crystals, characteristic of soft matter physics, this 
is no longer true because a number of simple constituents can share the same site. In real 
systems through surface adjustments and the spontaneous occurrence of defects it is possible 
to reach the equilibrium minimum free energy, that optimizes the ratio N/N s between the 
number of lattice sites N s , and the number of particles N. Conversely in simulations, 
periodic boundary conditions impose a constraint on the number of lattice sites N s , and 
at fixed number of molecules N, volume V and temperature T this prevents to reach the 
equilibrium state with the corresponding equilibrium ratio N/N s . In thermodynamic terms, 
the free energy cost to add or remove a lattice site, that we denote with /j, c , is too high to 
observe adjustments on typical simulation time scale. The free energy of the constrained 
system with fixed N s is determined not only by the conjugate couples (S, T), (II, V), (//, N), 
where S is the entropy, II the osmotic pressure and \x the chemical potential, but also by 
the site chemical potential /i c and its conjugate variable N s . In the constrained system, an 
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infinitesimal variation of F can be written as: 



dF(fi c ) = -S(fx c )dT - U(fi c )dV + n{iJ, c )dN + p c dN s (3) 

The expressions for the main thermodynamic constrained properties follow from the previous 
equation: 



dF 
dV 



N a NT 
Wiv ' N S VT 



The thermodynamic derivates, like the pressure and chemical potential, are performed keep- 
ing N s fixed, and this leads to consider constraint to the other properties. However 
N s , as well as fi c , is an unphysical variable and the results of the constrained system will be 
representative of the unconstrained situation only when site chemical potential /i c vanishes. 
Incidentally this condition restores the extensivity of the Gibbs free energy with the number 
of molecules in the system. 

G = F + PV = nN + fx c N s -> n c = (5) 

This condition is equivalent to say that the free energy cost to add or remove a site at 
equilibrium is zero and corresponds to the equilibrium value of the constraint. From a 
computational point of view the locus of points /i c (p, T) = 0, is obtained through canonical 
simulations at fixed densities p and temperatures T, varying both iV and V to find the 
specific values N, V for which: 

f v l\ s 

In order to determine the values (N, V) three independent measurements of F, II, [i are 
necessary. Following this formalism, we can recover predictions on the equilibrium uncon- 
strained system from a constrained simulation imposing eq. (|6|): 



U(N, V, T) = n(// c = 0) = U(N, V, T, N!«(N, V, T)) 
»{N, V, T) = fx(fx c = 0) = fi(N, V, T, N!«(N, V, T)) (7) 
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IV. COMPUTATIONAL ASPECTS 



For each thermodynamic state (p, T) it is necessary to perform several simulations in 
the constrained ensemble N,V,N S ,T, varying N/N s and correspondingly V/N s and com- 
puting F, /i, II until eq.([6]) is verified within the computational uncertainties. The estimate 
of the free energy, the osmotic pressure and the chemical potential, is achieved through a 
thermodynamic integration, the virial expression and the Widom insertion/removal method 
respectively |37J. The possibility to use the Widom insertion method in a dense phase is 
guaranteed by the softness of the interaction potentials. The implementation of Kirkwood 
coupling constant integration method [39] is significantly different from the usual Frenkel 
Ladd met hod [3^ ■ m the micellar crystals polymers can migrate from one aggregate to the 
next and therefore a molecule cannot be assigned univocally to a specific lattice site. To 
preserve such a feature along the coupling constant integration path, it is convenient to as- 
sume as reference system an ideal gas of dumbbells in an external field Q ex t that mimics the 
desired lattice geometry. The form of the external potential must allow the single molecule 
to jump from one Wigner-Size cell to the next as it is observed in the physical system. 
This can be achieved by summing over lattice sites, a dumbbell-site central pair interaction 
which vanishes at the edge of the WS cell. To compute the free energy difference between 
the original and the reference system, we employ a linear parametrization: 



where \J intra = Y^i Tab^) is the intradumbbell potential energy, U int er = J2i<j{ v AB+VBB)ij 
is the interdumbell potential energy and Q ex t = $a + $b is the energy of the external field 
acting on the A and the B blocks. The free energy is then obtained as |39j: 




where A = 1 corresponds to the reference system. The dumbbell-site interaction is in princi- 
ple arbitrary but for reason of efficiency it must be chosen to better reproduce some property 
of the original system. In the present case we have taken external fields acting on both A 
and the B blocks in such a way to roughly reproduce the density profiles, pa{t), Pb{^) within 
a micelle as obtained in a crystalline state simulation of the orginal system of interacting 
dumbbells. Strictly speaking, the inversion procedure from density profiles to single-particle 



U(X) = U intra + (1 - \)U inter + A$ 



ext 



(8) 
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FIG. 2: Radial density of core pa and corona pb (left hand), and effective external field (f>Aj<t>B 
(right hand) as function of the reduced distance R/Rg. a),c) density profiles for the real and the 
effective system for p/p* = 4; b),d) corresponding external fields, the vertical dashed line represent 
the potential cutoff r c = 2.3. 

interaction, requires the solution of two coupled integral equations due to the intradumbbell 
correlation. We have used a simplified approach based on simple barometric relations: 



Mr a ) = -K 3 log(p a (r a )) + Q a (10) 

where a = A, B and Q Q and X a are additive and dimensional constants respectively. In 
order to form a smooth function through the three dimensional space, it is important that 
each potential and its first derivative vanish at the edge of the Wigner-Size cell r c . To this 
aim we regularize and cut the above forms as 

ip a (r) + ip a (2r c - r) - 2tp a (r c ) r < r c 
r > r c 

These potentials, reported in fig. [2j represent the effective interactions between one lattice 
site and a block of a generic dumbbell. The external field to be used in our coupling 
constant integration is built by summing all contributions from different nodes of a lattice 
of the chosen symmetry {R m }, m = 1, N B : 

N s 

® a {r a ) = ^<p a {r a - R m ) a = A,B (12) 

m=l 
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The free energy of the reference system, i.e. a system of non interacting dumbbells in 



the lattice field of eq.(12) cannot be computed analitycally, but can be recast in the form 
of a statistical average and computed by MC simulation, as detailed in Appendix I. The 
numerical estimate of the integral over A in eq. ^ is obtained through the Gauss-Legendre 
quadrature scheme as detailed in Appendix II. 



V. DISPERSION OF MICELLES AND OD TRANSITION 

In the present work we limit our investigation to a single value of / = 0.6. The main reason 
is the large computational demand of the present technique. The other, more fundamental, 
reason is the limitation of the present dumbbell model, based on the zero density potential 
approximation, to investigate the intermediate density regime. Therefore this work is mainly 
meant as a methodological development to prepare the stage for studying the more accurate 
multiblob representation of diblock copolymer solutions [20J. 

We have performed MC simulations within the canonical ensemble, at various densities 
in the range p e [3.0,5.0]. In a previous work (23] the density of the onset of ordering was 
observed to decrease with increasing f, and we have chosen f=0.6 in order to explore the 
solid region without stretching too much the transferability hypothesis on the potentials. 
Clusters structure is characterized through the study of the distribution function P(n) that 
describes the probability of finding a cluster of size n. Above the CMC, this quantity shows 
a clear bimodal behavior with an exponential peak at the origin, signature of dispersed 
dumbbells, and a gaussian-like peak at a finite values of n signaling the presence of well 
defined clusters (see Figj3]) . The width of the cluster peak decreases with increasing density, 
meaning a smaller polydispersity and more stable aggregates at higher density. Note that 
P(n) has a minimum at n ~ 20 so that we consider an aggregate to be a cluster only when 
it includes more than 20 molecules. In Fig. [4] we report the average cluster size, computed 
as the ratio between the total number of molecules belonging to clusters (N — N na ) and the 
average number of detected aggregates (N m ) (panel a), and the fraction of non-aggregated 
dumbbells (N na )/N (panel b) as a function of the reduced density. Note that results at 
different densities are obtained with different system sizes and the smoothness of the behavior 
confirms the uniqueness of the equilibrium state for the system in the disordered phase, in 
contrast to the case of the ordered state (see below). We observe a roughly linear increase of 
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FIG. 3: Probability of a cluster of n molecules P(n), for different densities above the CMC. 
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FIG. 4: Panel a) estimate of the mean aggregation number of micelles as function of the reduced 
density; panel b) number of dispersed dumbbells as function of the reduced density on a semi-log 
scale, the dash-dotted line represents an exponential fit to the data. The dashed vertical lines 
indicate the estimated freezing density. 

the cluster size (N — N na ) / (N m ) up to the value p/p* = 3.8 above which the slope decreases 
considerably. This change of behavior signals the freezing transition for the micelles. The 
same indication can be obtained from the structure factor between the centres of mass of 
micelles 



S(k) = ( 



1 



N, 



cluster 



E 

^3 



ik-(R,* m Rem) 



(13) 



reported in fig. [5j where we see that at pj p* = 3.8 the amplitude of the first peak exceeds the 
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FIG. 5: Micelle-micelle structure factors S(k) as function of the reduced density above the CMC. 
For p/p* = 3.8, the first peak exceeds the Hansen- Verlet threshold (dashed line). 

Hansen- Verlet threshold value for freezing [38]. This criterion, proposed for simple fluids, 
it is found to hold also in the present case of micellar fluids. Indeed beyond p/p* = 3.8 we 
observe a slowing down of the micellar motion due to the sampling, despite the fact that the 
efficiency of the Monte Carlo sampling of single molecular moves is not reduced. Panel (b) 
of fig. [4] reports an exponential decrease of the number of non-aggregated molecules with 
density with only a marginal deviation above freezing. 

To further characterize the structure of micelles we report in fig. |6](a) the average densities 
of the micellar core, corona and of the dispersed particles as a function of the reduced density 
of the system. To compute these densities we have estimated the average volumes of the 
aggregates, both core and corona, by their radii of gyration. While the corona density 
remains always ~ 1 (in units of number of molecules over Rg) the core density is already 
~ 4 at the CMC and increases up to ~ 8 at p/p* =4. This is an anomaly of our present 
model due to the absence of any repulsions between A-blocks. The increase of core density 
up to p/p* ~ 3.5 is mainly due to the linear increase of the aggregation number (see figSJ^a)) 
not compensated by the almost linear increase of the core radius of gyration reported in fig. 
[6]^b). Above p/ p* ~ 3.5, the micellar core radius is observed to contract upon increasing the 
overall density, providing an enhanced increase of core density, even beyond the freezing point 
above which the aggregation number is almost constant (see figEfa)). The same qualitative 
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FIG. 6: Panel a) Densities of the core pa the corona ps (left axis), and of the dispersed dumbbells 
p na (right axis) as function of the reduced density; Panel b) Micellar Core and Corona radius as 
function of the reduced density. For sake of clarity the core radius is shifted by +1. The dashed 
line represent the estimated freezing density. 

behavior is observed for the radius of gyration of the corona although the turning point is 
at p/p* ~ 3.7 and the following decrease is less important. Finally in fig^a) we also report 
the density of the dispersed molecules computed as the number of dispersed molecules over 
the volume available to the solvent. The latter has been obtained using the corona radius of 
gyration to estimate the volume occupied by the aggregates. A linear decrease is observed, 
meaning that the decrease of available volume partially compensates for the exponential 
decrease of available molecules observed in fig[4^b). 

In contrast to low density disordered states, the "experimental" observation is that be- 
yond the freezing point the stationary state depends considerably on the total number of 



molecules in the system. In particular for N in specific ranges (see section VI) the system 
is observed to spontaneously order in defective crystalline structures of micelles. As an 
example, at p/p* = 4 a reasonable guess for the number of dumbbells iV compatible with 
the number of micellar sites for a BCC cluster lattice, allows us to observe the spontaneous 
emergence of a pattern of peaks in the static structure factor reminiscent of a BCC order- 
ing. Conversely, for an arbitrary number of dumbbells the system shows a structured S(k) 
with small peaks and slow structural adjustments, signature of an arrested glassy phase. 
Simulations beyond the estimated freezing density for different numbers of molecules reveal 
a significant change of the P(n) cluster peak, associated with the occurrence of a sponta- 
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FIG. 7: Comparisons between the probabilities of finding a cluster of size n, P{n) (left panel), and 
between the micellar radial distribution functions g{r) (right panel) at the same density p/p* = 4. 
and different number of dumbbells N = 4500, 2500. 

neously ordered phase. An example is given in figure [7] where we compare P{n) at p/p* = 4 
for a BCC phase and a disordered phase. A decrease in the position of the peak and in 
its width in the ordered phase are evident, accompanied by the appearance of a splitting of 
the second peak in the radial distribution function between the centres of mass of micelles. 
This phenomenon is associated with a reduction of the system polydispersity as counterpart 
of the change of local and long range order in the crystalline state. This is an aspect of 
the complex interplay between the long range ordering of the aggregates and the micellar 
internal structure. 



VI. FREE ENERGY CALCULATIONS 



In this section we apply the thermodynamic formalism developed previously to the free 
energy calculations of cluster crystalline structures of different symmetry such as simple cubic 
(SC), BCC, FCC, A15 and diamond, and at two densities p/p* = 4,5. The preparation of 
the system in specific crystalline structures has been achieved by applying the same external 



field described in section [TV] for the A integration procedure but keeping at the same time the 
intermolecular interactions active. It turns out to be possible for a given density pj p*, crystal 
symmetry and number of lattice sites N s , to find a window for the number of dumbbells 
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N in the system, in which the different clustered phases are dynamically stable during the 
MC simulation. This is true in all cases but the diamond structure for which, even with 
the field turned on, the clusters were not found to order according to the symmetry of the 
external field. The reason is that in the diamond structure the nearest neighbors distance is 
shorter than in the other structures and would impose the formation of micelles too small 
to be stable. Even with the external field of reasonable amplitude turned on, the system 
prefers to occupy only one site over two and form larger aggregates. Upon turning off the 
field, the SC structure shows a dynamical instability with a rapid melting independently 
of N. The other three structures were dynamically stable and for each of them we have 
performed several simulation with different (N,V) in the constrained ensemble (p/p*,N s ), 
namely at fixed density, lattice geometry and number of lattice sites. In order to compute 
the free energy of each thermodynamic constrained state (N, V, N s ) by the coupling constant 
integration procedure, see eq.([9]), further simulations are necessary for different values of the 
coupling parameter A chosen according to the Gauss-Legendre method. We have used 8 
quadrature points and we report the numerical details in appendix I. 

In tables [I] and [TT] we collect the results for the excess free energy, osmotic pressure and 
chemical potentials, including p c , for different structures at the densities p/p* = 4 and 5, 
respectively. The first observation are the large uncertainties on the determination of p c . 
This is due mainly to the large size of our aggregates since, according to eq. ([6]), p c is obtained 
as combination of specific properties, whose uncertainties are quite small, multiplied by the 
average number of molecules per lattice site N/N s . It is then clear that uncertainties on p c 
are two orders of magnitude larger than uncertanties on any other intensive thermodynamic 
property. This fact makes is quite difficult to obtain an accurate determination of the number 
of particles at which p c vanishes. 

At p/p* = 4, excess free energies and chemical potentials are largely insensitive to N/N s 
and to the lattice geometry. Conversely the excess osmotic pressure decreases for increas- 
ing N/N s probably due to its dependence on the lattice spacing which also increases with 
N. Despite the large uncertainties on p c , our results allow to identify the free energies of 
the unconstrained crystals by a simple linear interpolation between positive and negative 
values of p c (for the A15 structure we assume the point at p c = 0.07 ± 0.9 to be repre- 
sentative of the unconstrained system). The interpolated values for the excess free energies 
per dumbbell, which allows to determine the most stable structure at the given density, 
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are fncc = 7.6767(5), /ais = 7.6776(6) and Jfcc = 7.6792(4). Therefore our procedure 
indicates the BCC structure to be the most stable among the three considered, with the 
A15 structure only slightly above and within the error bars. These are both open structures 
in contrast to the close packed FCC structure which appears to have an higher free energy. 
Note that the variation of the Helmholtz free energy among the various structures is of the 
same order of magnitude as the observed variations with N/N s within the same structure, 
so that it would have been impossible to determine the ordering of the structures without 
determining the unconstrained state at /i c = 0. 

At p/p* = 5 a similar study for BCC and FCC lattice, as summarized in table |H| , 
shows a clear dependence of the excess free energy, as well of the osmotic pressure and 
fi c , on the number of molecules. However as at lower density, the excess chemical potential 
is almost insensitive to the number of particles. Even in this case the uncertainties on 
fi c are large but the stability of the crystalline structure for systems with both positive 
and negative values of [i c allows easily to identify the equilibrium unconstrained state by 
linear interpolation. The extracted values for the excess free energies per dumbbell are 
Jbcc — 9.12(1), Jfcc — 9.184(3) indicating again the BCC structure as more stable for this 
system. 




n 



FIG. 8: BCC structure: comparison of the cluster peaks of P(n) for p/p* =4 (left) and p/p* = 5 
(right) at different numbers of dumbbells. The dashed lines represent gaussian fits to the peaks. 

The present procedure to determine the unconstrained equilibrium state, allows to discuss 
the true equilibrium properties of the micelles in the symmetry broken phase. For the stable 
BCC phase, the variation of P(n) with the total number of molecules is reported in figure 
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[8] for the two densities investigated. At both densities a shift of the cluster peak to higher 
aggregation numbers for increasing N is observed. The shape of the peak is fairly well 
fitted by a gaussian function with variance independent of N. The observed decrease of the 
amplitude of the peak with N is balanced by an observed increases of the number of dispersed 
molecules. Numerical results for the mean value and the variance of the fitted gaussian are 
reported as n c and its uncertainty in tables [l| and |H| for p/p* = 4 and 5, respectively. The 
interpolated equilibrium values for the BCC structure are n c = 81(12),n c = 115(12) for 
p/p* =4 and p/p* = 5 respectively. 

Tables [II and [TT] also report the values of the Lindemann ratio L and of the average number 
of vacancies per lattice site, (N v )/N a , for the various structures as a function of N. In general 
for a given structure, both quantities decrease for increasing N meaning that the mechanical 
stability of the lattice increases with N. However this behavior is limited to a finite range 
of N values above which the crystal structure is found to spontaneously melt. At p/p* =4 
the Lindemann ratio is significantly larger than the classical value at melting (£ ~ 0.15) 
and found to be £/ cc ~ 0.25 and £f, cc ~ Cais ~ 0.29 at the unconstrained state. The 
equilibrium value 0.29, reminiscent of the value observed for delocalized quantum particles 
near melting|40], suggests that the internal micellar structure plays the role of quantum 
fluctuations with respect to the vibrations of the centre of the composite particle. Note that 
we have not computed the free energy of the disordered liquid state at p/p* =4 but the 
BCC crystal stability against the liquid phase is supported by the observed spontaneous 
crystallization of a sample with N = 4500. According to table [I] this system has \i c ~ 1, 
resulting in a free energy higher than the unconstrained state, and therefore proving the 
stability of the latter with respect to the disordered phase. At p/p* = 5 the same qualitative 
behavior is observed with smaller values of £'s and (N v ) /N s for all structures meaning that 
this point is deeper inside the crystal phase. At this density and for the stable BCC phase, 
we also report in figure [9] the peaks of the micellar S(k) which allows to directly extract 
the Debye- Waller factor W j4T]. We observe a trend of W compatible with the increased 
stability with N. Note that the equilibrium state does not corresponds to the mechanically 
most stable state. This is natural since maximal mechanical stability corresponds to the 
state with no vacancies while the thermodynamic equilibrium state has a non vanishing 
concentration of vacancies. 
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FIG. 9: Logarithm of the micelle-micelle structure factor S(k) versus the square modulus of the 
reciprocal vector multiplied by the lattice spacing a, for a BCC structure of clusters at p/p* = 5., 
for an increasing number of molecules. The dashed line indicates the onset of anharmonicity. 

TABLE I: Main thermodynamic and structural properties for different structures and number of 
dumbbells for p/p* = 4. and / = 0.6. N s , n c , C and {N v )/N s represent the number of sites, 
the mean aggregation number of micelles, the Lindemann ratio and the equilibrium fraction of 
vacancies respectively. 



¥ S -Structure 


N 


N/N s 


n c 


F ex /N p ex 


U ex /p 


Pc 


Ap c 
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2916 


91 


85(12) 


7.6794(4) 13.99(1) 


6.3031(8) 


-0.7 
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6.382(1) 
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0.260(2) 0.0422(9) 
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7.6776(4) 14.00(1) 


6.323(1) 


0.07 
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0.290(2) 0.061(1) 
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TABLE II: Main thermodynamic and structural properties for different structures and number of 
dumbbells for p/p* = 5. and / = 0.6. The definition of the symbols is the same as in Table [I] 
^-Structure N N/N s n c F ex /N p ex U ex /p p c Ap c C (N v )/N s 
32 -FCC 4300 134 131(11) 9.1881(8) 16.21(1) 6.9836(6) -5.1 1.3 0.06889(2) 0.0003(1) 
32 - FCC 4000 125 122(11) 9.1813(8) 16.20(2) 7.0477(6) 3.6 2.5 0.0747(8) 0.0008(1) 
32 - FCC 3680 115 112(11) 9.1667(8) 16.26(2) 7.1320(4) 4.4 2.3 0.264(3) 0.060(1) 
32 - FCC 3500 109 107(11) 9.1703(8) 16.29(2) 7.1842(4) 7.0 1.6 0.291(2) 0.062(1) 
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120 
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111 


109(11) 
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BCC 
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101 


100(11) 


54- 
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94 
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9.0885(9) 16.21(2) 7.0722(4) 
9.1807(9) 16.25(2) 7.1026(4) 
9.1895(9) 16.28(2) 7.2273(4) 
9.1985(9) 16.38(2) 7.3566(4) 



-5.9 2.4 0.083(2) 0.0017(2) 

13.7 2.2 0.0901) 0.0020(2) 

13.9 2.0 0.103(1) 0.0037(2) 

16.5 1.9 0.124(3) 0.0070(5) 



VII. CONCLUSION 



We have reported a quantitative investigation of the relative stability of self- assembling 
miceliar crystalline phases of diblock copolymer in semi-dilute solutions. We have considered 
an athermal ISS model in which a copolymer has a block (A) in theta solvent conditions 
while the other block (B) in good solvent conditions, and the A-B interaction are mutually- 
avoiding. We have adopted the minimal coarse-grained model for such systems, consisting 
of mapping a long diblock copolymer onto a soft-core dumbbell. The intra- and inter- 
molecular interactions are obtained in the infinite dilution limit and then used at finite 
dilution on the basis of an untested transferability assumption. Beyond a CMC, the system 
of dumbbells is found to self-assemble into spherical clusters which form a liquid-like system 
of micelles. Upon further increase of density a ODT between the miceliar liquid and a 
miceliar crystal is met. Our investigation has aimed to determine the most stable crystalline 
state among several candidate cubic structures. This program requires the calculation of 
the free energies for the various candidate structures and has been performed by applying 
to our system a recently proposed method to compute free energies of multiple-occupancy 
crystals [2"4ll2"o] . Our study has been limited to a single value of the relative block length ratio 
/ = Ma/ {Ma + Mb) = 0.6 and to only two densities p/p* = 4 and 5, the former close to the 
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(still undetermined) transition density, the latter more inside the crystal phase. The reason 
is partially the large computational cost of the present method, but more fundamentally the 
qualitative character of our minimal coarse-grained model to represent diblock copolymer 
solutions at intermediate density. Therefore the present study is mainly methodological and 
has prepared the stage for the application of the free energy technique to the more accurate 
multi-blob representation [19"} I2"U]. 

We have performed MC simulations for increasing density in the range p/p* G [3,5] 
with the aims of (i) characterizing the micellar structure with density in the liquid phase, 
an aspect which has not been fully reported in previous works (221 123] , and (ii) finding the 
freezing density at pj p* f» 3.8, as signaled by the Hansen- Verlet criterion for the S(k) and by 
a marked change of behavior in the mean aggregation number. Furthermore at p/p* = 4 and 
5, for several crystalline symmetries, we have performed a number of canonical simulations 
for increasing N/N s in order to determine the equilibrium state, at which the site chemical 
potential p c vanishes. This is the unique equilibrium symmetry broken state of the system 
that needs to be known in order to compare different crystal symmetries without any bias 
imposed by the simulation constraints. We have found that our present system prefers the 
less packed BCC structure over the more packed FCC ones, with the A15 structure, an 
open non-Bravais structure which minimize the surface to volume ratio of the Wigner-Size 
cell [2H1 El] , almost degenerate with the BCC one. At p/p* = 5, well inside the crystal phase, 
the ratio of the corona layer thickness to the core radius L / R c , both estimate from the radii 
of gyration of core and corona, is L/R c = 1.3(4) in good agreement with the experimental 
value L/R c = 1.2 for poly(stirene-b-isoprene)22-12 with f=0.6 (see table I of |6j). Also 
in qualitative agreement with the experiments, we find that our system prefers the BCC 
symmetry over the FCC one. Knowing the true equilibrium state of the system we are in 
the position to compute structural equilibrium properties such as the Lindemann ratio C, 
the average number of vacancies per lattice site (N v )/N s and the Debye- Waller factor W. 
Near melting we observe a large Lindemann ratio (~ 0.29) reminiscent of the value for low 
temperature quantum particles f4U], and a substantial number of vacancies per lattice site 
(~ 6%). Both quantities decrease when going deeper inside the crystal phase and we observe 
a residual number of vacancies per lattice site of ~ 0.2% at p/p* = 5. 

One limitation of our present approach is related to the "small" size of the systems 
considered. Indeed we examine systems of small number of micelles (32, 54 and 64) and we 
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can expect a finite size effect on the computed free energies. Such effects wili not change the 
ordering of structures far from the melting transition but could be relevant for an accurate 
location of the melting line and, because of the small free energy differences observed near 
melting, to determine the stable structure upon crystallization. A possible route to address 
this problem would be to invoke a second level of coarse graining in which each micelle is 
replaced by a point particle and the density dependent pair potential between two micelles is 
extracted by the radial distribution function between the centre of mass of micelles computed 
in the present dumbbells coarse-grained representation [221 EZ]- With this further coarse 
graining, roughly 100 dumbbells are mapped on to a single particle with an important gain in 
efficiency and the possibility to study systems of thousands of micelles with negligible finite 
size effects. Despite a number of approximations involved in this further coarse-graining 
step, namely the neglect of the polydispersity of micelles and of the exchange of polymers 
between micelles and with the dispersed phase, this procedure was found to reproduce very 
accurately the structure of the micelles in the liquid phase and its extension to determine 
the relative stability of several crystal phases has been attempted [27J. It was found a general 
preference for the A15 structure followed by the FCC structure, while the BCC structure was 
found to be only marginally stable. However the micelle-micelle potential can be extracted 
from the structure, only in the continuous symmetry phase and its determination at solid 
densities was limited by the possibility to prevent the system from crystallizing. We have 
seen here however that the structure of micelles changes significantly upon crystallization at 
fixed density, which would imply that the effective micelle-micelle potential would depend 
not only on the density but also on the symmetry of the state. Unfortunately a procedure to 
extract a pair interaction from structures with broken symmetries has still to be devised and 
would require the extension of integral equation theory to the crystalline phase. We believe 
this limitation in obtaining the pair potentials is the main reason for the disagreement with 
our present results. 

Very recently the extension of the multi-blob approach [19] to diblock copolymers has 
been proposed and applied to study the phase diagram of the ISS model for various values 
of / and at various densities [20] • In the crystalline phases, the same dynamical stability in a 
limited range of N/N s as observed by us, has been reported and the Helmholtz free energy 
per polymer was found to be independent on N/N s . On the basis of this independence, the 
authors concluded that the system is not able to select an optimal cluster size and therefore 
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an optimal lattice spacing. However no evaluation of the unconstrained equilibrium state 
at He — has been attempted. From our present results we see that F/N and fi are fairly 
insensible to N/N s . However this is not proving the inability of the system to select a 
specific aggregation number and lattice spacing, but it is signaling the presence of the extra 
constraint imposed by the fixed number of lattice sites in the simulation box. We believe 
our present analysis, when applied to the multiblob representation, will locate univocally 
the unconstrained equlibrium state and its properties. 

It is well know that temperature is a relevant variable in determining the phase diagram 
of experimental diblock copolymer solutions. Changing temperature the selectivity of the 
solvent to the two blocks of the copolymer can be continuously varied, giving rise to ther- 
motropic transitions and to a rich phase diagram. Our present model, being athermal, is 
not capable to reproduce this phenomenology. The extension of the present coarse-graining 
strategy to the case of homopolymers in solvent of varying quality has already appeared[42j. 
At variance with the case of chains in good solvent, it was observed that the density depen- 
dence of effective potentials between polymers is quite significant and becomes increasingly 
important when approaching the theta point. This would imply that a minimal dumbbell 
model of diblock copolymer solutions with zero density effective potentials is a poor repre- 
sentation of the system. However, the temperature effect could be included in the present 
coarse-graining strategy by adopting the multiblob representation based on zero density po- 
tentials and varying the level of coarse-graining according to the chosen temperature. Work 
along this direction is in progress. 

Finally we briefly discuss the relation between the present coarse-graining strategy and the 
Self-Consisted Mean Field Theory (SCMFT) based methods [TTHT5] . The present simulation 
approach, although based on coarse-grained models, remains computationally much heavier 
than SCMFT. Moreover determining free energies requires special methods as shown in the 
present work, and inferring accurate phase boundaries requires a method to estimate finite 
size effects still to be developed. Therefore at present SCMFT methods are superior in 
their ability to reproduce the entire phase diagram with limited effort. Nevertheless, our 
multiscale strategy to go from the full monomer system to the system of many micelles 
in a consistent and controlled way provides detailed information on new phenomena such 
as the interplay between the internal structure of the micelles and their order which has 
not been addressed previously. This interplay can determine a dependence of the effective 
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micelle-micelle interaction on the ordering of the system which can modify the predictions 
of the phase diagram based on the effective micelle represent at ion [T2j 127] . 
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VIII. APPENDIX I 



In this appendix we describe the procedure to compute the free energy of the reference 



system, i.e. a system of non interacting dumbbells in the lattice field provided by eq.(12). 
The partition function of this reference system can be expressed as an ensemble average. 
To this end it is convenient to write the single molecule hamiltonian in terms of molecular 
coordinates and momenta (r, R), (p, P) 

N 



Pl + F f + iA B {n) + $ A (ifc - |) + Sato + 1; 



(14) 



where (r,p),(R,P) are the internal and center of mass position and momentum of the 
dumbbell defined respectively as: 

D r A + r B 

R= —r~ 

r =r A -r B 

P =2mR 
mr 

P ( 15 ) 

with ta and r B being the coordinates of the A and B blocks of the generic dumbbells 
respectively. From now on we set m = 1. Tracing out the momenta the partition function 
reduces to the product of the ideal gas term and an average over the equilibrium bond 
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distribution: 



• Qn(V,T) 



1 ? dr N e -PY.f[lAB(n)+^A(R 1 ~r t /2)+^ B {R 1 +r l /2)} 



A^AfiV! 



n N f r p-MA B {r) \ ^ 

( / ,/*?,/»._ -/3[* i4 (fl-r/2)+* B (fl+r/2)] 

A^AfiVlU q mtra 

N / r \ N 

Qintra / / J L> / p -P(&a+&b)\ . 



dfl(e-«^ + *^) 1WrB (16) 



where qmtra — J dre~^ AB ^ is the configurational partition function of a single dumbbell, 
Acm and A M are the de Broglie thermal wave lengths of the center of mass and of the internal 
coordinate respectively, and (. . . ) intra indicate an average over the bond distribution funtion 
Pintra( r ) = e~ '^ AB ^ j 1 q intra . The reference free energy follows as: 



F(X = 1)/N = -k B T\og [jrj^ ~ k B T\og^- k B T\og (I J dR^^) Pmtr )j 
= F ideal (N,V,T)/N-k B Tlog(^ J dR(e-^ + ^) Pi A (17) 

To obtain the excess free energy it is necessary to perform a MC simulation of a single 
dumbbell in a volume V at temperature T, sampling the bond vector according to Pmtraif) 
and averaging along the stochastic trajectory the quantity e -/ 3 [^A(R-r/2)+<s> B (R+r/2)} _ p Qr a 
faster convergence center of mass moves with amplitude Va /2 and unitary acceptance has 
been implemented. 



IX. APPENDIX II 



In this appendix we report details of the application of the Gauss-Legendre method to 
compute the A-integral in eq. ([91). At each A value a N, V,T simulation is necessary to 
compute {§ ex t — Ui n ter)\- The softness of the potentials allows us the use of standard single 
particle moves with large amplitude and hence a good sampling of the previous expression at 
finite A. This is no longer true when A — > 0, because the vanishing of the pinning potential 
in the Boltzmann weight, progressively allows a slow drift of the system with respect to the 
location of the field resulting in large fluctuations and slow convergence of the term (& e xt)\- 
To improve the statistics at A = we have randomly sampled the geometrical center of 
the field inside the simulation box. For small but finite A-s (~ 0(1/N)) this procedure is 
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• Gauss- Lobatto interpolant 
□ Gauss-Legendre interpolant - 
— 3rd degree correction polynomial 



FIG. 10: Thermodynamic integrand as function of the coupling parameter A for a BCC structure 
and p/p* = 5. Closed circles indicates the Gauss-Lobatto points, while the open squares the Gauss- 
Legendre points. The solid, dashed and dot-dashed lines correspond to the 3rd order correction at 
small A, to the Gauss-Legendre and to the Gauss-Lobatto interpolants respectively. 

doomed to fail and a Metropolis sampling of the centre of the field would be necessary but 
much slower. In figure 10 we show a typical behavior of the integrand in eq.([9]). We note a 



very smooth behavior with A except at small A where a sharp increase is observed due to the 
progressive vanishing of the pinning potential. To accurately integrate this behaviour with 
a reasonable effort we can either apply the Gauss-Legendre scheme or the Gauss-Lobatto 
scheme [13], the difference being that in the former scheme only interior points to the inter- 
val are used, while the latter scheme employes also the extremes of the interval. In order 
to compare the accuracy of the Gauss-Legendre scheme and the Gauss-Lobatto one with 8 



and 10 points respectively, we report in figure [10] the associated interpolating polynomials. 
We note that the interpolating Gauss-Legendre polynomial well represent the expected be- 
haviour, except in a small interval below the first point Ai where we expect a sharp increase 
toward the value at A = to occur. Conversely the Gauss-Lobatto interpolating polynomial 
exhibits an oscillating behaviour which is unphysical. For this reason we have adopted the 
Gauss-Legendre scheme and implemented a modification to approximately correct for the 
missing area near the origin. We model the unknown integrand function between the origin 
and the point Ai by a third order polynomial uniquely defined by matching the values if 
the integrand function and its first derivative at A = and at A = 1. The value of the first 
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derivative at the origin can be estimated by a simple cumulant expansion of the integrand 

(AU) X = ($ ext - U mter ) x ~ (AU) - A ([AC/ - (AU) }\ (18) 

while the derivative in Ai can be obtained as the derivative of the Gauss-Legendre interpolat- 
ing polynomial. Within this scheme, the corrected estimate of the thermodynamic integral 
is 

8 

AF = J2(^U}xM^) - WO, Ai) + WO, Ax) (19) 

where the first term is the Gauss-Legendre quadrature formula, the second term is the 
numerical integral of the Gauss-Legendre interpolating polynomial between [0, Ai] and the 
last term is the integral over the same interval of the new interpolating polynomial. This 
correction, of order Ai(AC/)o, is essential to obtain an unbiased estimate of the free energy 
and in our present case turns out to be roughly 0.1% of the value of the total free energy, 
an order of magnitude larger than the statistical error. 
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